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ABSTRACT. We study the probability of rain before time t for the family of tempered 
stable Levy insurance risk processes, which includes the spectrally positive inverse Gauss- 
ian processes. Numerical approximations of the rain time distribution are derived via the 
Laplace transform of the asymptotic rain time distribution, for which we have an explicit 
expression. These are benchmarked against simulations based on importance sampling 
using stable processes. Theoretical consequences of the asymptotic formulae are found 
to indicate some potential drawbacks to the use of the inverse Gaussian process as a risk 
reserve process. We offer as alternatives natural generalizations which fall within the tem- 
pered stable family of processes. 



1. Introduction 

The risk reserve of an insurance company has traditionally been modelled as a com- 
pound Poisson process with drift. In recent years more general Levy processes have been 
proposed, among them the inverse Gaussian family of processes. Such processes have 
been found to approximate reasonably well a wide range of aggregate claims distributions 
[9]. While the probability of eventual ruin has received a lot of attention, arguably of equal 
importance in practice is the probability of ruin before some finite time horizon. Our paper 
aims to study the probability of ruin before time t for the inverse Gaussian family and a 
natural generalisation, the tempered stable processes. 

The basis of our investigation is the recent asymptotic representation, as the initial re- 
serve grows large, of the ruin time distribution for more general "medium-heavy" convo- 
lution equivalent Levy processes [19, 21]. This representation, via the calculation of its 
Laplace transform, lends itself to a numerical approximation of the ruin time distribution, 
which is then benchmarked against values obtained by simulation. Thus we are able to 
illustrate the use of a broad, relatively simple and computationally tractable family of pro- 
cesses with which to model the risk reserve process. 

We find that the asymptotic representation performs well even when the initial capital 
is relatively small, contrary to a view that asymptotic formulas may only be useful when 
the initial capital becomes extremely large. Additionally, the asymptotic representation 
provides some interesting insight with regard to safety loading management. When a real- 
istic safety loading is specified in the insurance risk model, we show that processes within 
the tempered stable family may exhibit undesirable exponential growth (in time) of the 
ruin probabilities, at least asymptotically. This indicates that some caution may need to 
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be exercised in the choice of model and to aid with this task, we derive a useful relation- 
ship between the parameters to avoid an unpleasant scenario. This might have interesting 
implications for practitioners concerned with safety loading management. 

Empirically we also observe that the asymptotic formula provides a useful lower bound 
for the ruin probability that can be combined with the infinite horizon ruin probability to 
provide a practical approximation of the true ruin probability. 

1.1. Levy insurance risk model. Let X = {X, : t > 0}, Xq — 0, be a Levy process defined 
on (Cl,^,P), with canonical triplet (yx, a^,XIx)- The characteristic function of X then has 
the Levy-Khintchine representation Ee' ex ' — e'^' 9 ', where 

^x(O) = iey x -\ole 2 + j (e iex -l-i0xl {H<1} )nx(ck), forQeM. (1) 

In the general Levy insurance risk model, the claim surplus process, which represents 
the excess in claims over income, is modelled by a Levy process X with X, — > — °° almost 
surely. Claims are represented by positive jumps, while premia and other income produce 
a downward drift in X. The insurance company starts with a positive reserve u, and ruin 
occurs if this level is exceeded by X. The assumption X t — » — °° a.s. is a reflection of 
the premium being set to avoid certain ruin. This setup generalises the classical Cramer- 
Lundberg model in which 

N, 

X t = Y,Ui-pt, (2) 

/=i 

where the nonnegative random variables f/, form an i.i.d. sequence with finite mean fx, N t 
is an independent rate X Poisson process, and p > Xji. Here t/, models the size of the /th 
claim and p represents the rate of premium inflow. The assumption p > Xjl is the net profit 
condition needed to ensure that X t — s- — °° a.s. See [2] for background. 

1.2. The convolution equivalent model. A natural class which includes the tempered sta- 
ble distribution and the inverse Gaussian distribution is the class of convolution equivalent 
distributions. Definitions and basic results for convolution equivalent distributions and the 
corresponding convolution equivalent Levy insurance risk processes are set out in detail 
in Kluppelberg, Kyprianou and Mailer [24] and Griffin and Mailer [21], and associated 
papers, so we only outline the main ideas here. A comparison of the medium heavy convo- 
lution equivalent condition, the light-tailed Cramer condition (Ee v ° Xi = 1 for some Vo > 0) 
and the heavy tailed subexponential condition can also be found in [21]. 

Denote the class of (non-negative) convolution equivalent distributions of index a > 
by ,y( a \ A Levy process is said to be convolution equivalent 1 , written 

X+ G ,y {a) for some a > 0, (3) 

if the distribution of is in J^W for some a > 0. The convolution equivalent Levy 
insurance risk model is one in which 

x + £ y(a) for some a > o and X, ^ -oo a.s. (4) 

Membership of S fi ^- a \ by definition, is a property of the positive tail of the distribution 
of X\ . Condition (3) can equivalently be expressed in terms of the positive tail ILt (u) = 
Tlx((u,°°)) of the Levy measure (see [24]). Assuming TI x (xq) > for some > 0, so 
that X has positive jumps with probability 1, we say that Tl x e ,y( a ) if the same is true of 



See Borovkov and Borovkov [6] and Foss, Korshunov and Zachary [16] for further background on subexpo- 
nential and convolution equivalent distributions. 
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the corresponding renormalised tail (Tl x (-)/Tl x (xq)) A 1. With this understanding, (3) is 
equivalent to 

€ J^ (a) for some a > 0. (5) 

Convolution equivalent distributions of index a have exponential moments of order a, 
but of no larger orders. Thus, if \j/x denotes the cumulant of X, so that 

Ee m =e tyx(P) j 

then Wx{P) is nn i te if an d only if /3 < a. 

Some asymptotic aspects of the model (2) where U\ has a convolution equivalent distri- 
bution were recently considered by Tang and Wei [33]. In particular, explicit asymptotic 
formulas for the Gerber-Shiu function in the infinite horizon case were derived. Theoretical 
and numerical comparisons between models under the Cramer condition or a convolution 
equivalent condition were recently carried out in [22] for general Levy insurance risk pro- 
cesses. It was observed that the "medium-heavy" regime transitions continuously into the 
"light-tailed" Cramer regime as certain parameters describing the models are varied. The 
convolution equivalent model was suggested as providing a broad and flexible apparatus 
for modelling the insurance risk process. 

1.3. Eventual ruin. Convolution equivalent Levy processes were introduced into risk the- 
ory in [24]. In addition to (5), [24] assumed 

Ee aXl < 1. (6) 

Condition (6) implies that {e aX ' ) t >Q is a non-negative supermartingale from which it fol- 
lows thatX r — > — °° a.s., so the second condition in (4) is automatic in this case. 
For a given initial reserve u > 0, the ruin time is defined by 

t(k) =inf{r >0:X t > u}. (7) 

The main results in [24] include the following asymptotic estimate for the probability of 
eventual ruin. Assume (5) and (6). Then 

15(h) -Yx(cc) 

where 

X, = sup X s . (9) 

0<s<t 

This expression for the limit differs in form from that given in [24], but is equivalent; see 
Remark 1. Under (6), yfx(a) < and Ee aY °" < oo. If Ee aXl e [1,°°) then Ee aX - = <», but 
Ee aX < < oo for all t > 0; see Lemma 2.1 in [19]. 

1.4. Ruin in finite time. A more difficult problem than the probability of eventual ruin is 
to find the distribution of the ruin time itself. For convolution equivalent processes, partial 
results in this direction were obtained by Braverman [7], Braverman and Samorodnitsky 
[8], and Albin and Sunden [l]. 2 More recently, the following explicitly defined asymptotic 
estimate was obtained in Griffin [19] and Griffin and Mailer [21] under the sole assumption 
(5): 

P{t(u) <t)^Ux(u)B(t)+o(Ux(u)) a.s. w^oo, (10) 



2 

Heavy tailed (subexponential processes) are treated in Asmussen and Kliippelberg [3]. For the light-tailed 
"Cramer case", see [2]. 
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where the function B(t) satisfies 

B(t)= I' ds. (11) 



o 

Under Condition (6), the estimate in (10) is uniform in t > 0, B is bounded, and by mono- 
tone convergence, B(t) increases as t —> °° to 

Ee a *°° 

BH = (12) 

which coincides with the limit in (8). In this case the estimate may be rewritten in a more 
intuitively appealing form. From (8), (10) and (12), we have for t > 

PMu)<t)=P(T(u) <<*>)( 2fL+o(l) ) asif^oo. (13) 

\B{oo) J 

Thus, asymptotically as u — > oo, P(x(u) < t) factors as the product of the probability of 
eventual ruin and a distribution function in t given by B{t)/B{°°). Moreover this estimate 
is uniform in t > 0. 

When Ee aX ' > 1, (10) provides an estimate which is uniform on compact sets in t > 0, 
but the function B is unbounded and the limit in (8) is infinite. Thus P(l(u) < t) is no 
longer proportional to P{x{u) < °°) and (13) does not hold. 

1.5. Overview. The quite explicitly defined form of B(t) in (11) opens the possibility 
of calculating it numerically for an appropriate class of models, with the hope that the 
estimates may be used for guidance in some real-life modelling situations so as to derive 
useful information about the ruin time distribution. But in order to implement this program, 
a number of questions must be addressed. First, what can be said about properties of Bl 
Second, how do we obtain a good numerical approximation for the expression given by 
(11)? Third, once numerical results are at our disposal, how well do the approximations 
(10) and (13) perform compared to, say, a direct simulation of the ruin time probabilities 
for different values of u and f? The aim of this paper is to give answers to these questions. 



2. Some Fluctuation Theory 

In order to investigate properties of B we need to introduce some notation and a few 
basic results from the fluctuation theory of Levy processes as set out in Bertoin [4], Sato 
[31] and Kyprianou [25]. 

2.1. Inverse local-time and ladder height processes. Let (L, r 1 ,H t ) t >o denote the bivari- 
ate ascending inverse local time and ladder height subordinator process of X. The bivariate 
descending inverse local time and ladder height subordinator is denoted by (L t: H t ) t >Q. 
Their Laplace exponents x{a,b) and K(a,b) are defined, for values of a,b £ R for which 
the expectations are finite, by 

e -K(",b) =E ( e -aLY 1 -bH l . l<L ^ and e -?M) =£ (,-^ 1 -M l;1< q i (U) 

The random variables Loo and Loo are exponentially distributed with parameters q > and 
q > respectively, with the understanding that if q or q is zero then the resulting random 
variable is identically infinite. We can write 

K(a,b) =q + d L - ia + d H b+ ( ( (\-e-"'- bi; )ll L - iH {&t,&h), (15) 

Jt>0 Jh>0 V ' 
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where d L -i > and d# > are drift constants, and Yl L -i H (dt ,dh) is the bivariate Levy 
measure of (L ,H). Similarly, 

K(a,b)=q + d ? ^a + dffb+ I [ ( 1 - e - a, - bh \ JJ ? ^ s(d/,d/i). 

Jt>0 Jh>0 V / ' 

These integrals are finite at least for a > 0, b > 0. We denote the marginal Levy measures 
of L~ l and H by II L -i (df ) and n#(d/z), and similarly for the corresponding hat quantities. 

When lim,^ooX r = — °° a.s. the increasing ladder process (L,~ ,H t ) t >o is defective. In 
that case (L~ l ,H) is obtained from a non-defective bivariate subordinator (J? -1 , Jf? ) by 
independent exponential killing with rate q > 0. The decreasing ladder process (L, . H t ),>o 
is proper when X r — > — °° a.s., and we then have = 0. 

The Wiener-Hopf factors of X may be expressed in terms of these Laplace exponents. 
In particular 

Ee~*> = ^ (16) 
K{o,a) 

where e is independent of X and exponentially distributed with parameter 8, and a > 0. If 
TT^ e y( a \ then Ee aXl < «>, £e affl < °o and (16) remains true for a > -a. 

The Wiener-Hopf factorization involves an arbitrary constant which we may take to be 
one by choice of normalization of the local times. In other words we assume 

-logEe' ex > = [-logEe' m ][-logEe- ieSi ]. (17) 



2.2. The spectrally positive case. When X is spectrally positive, that is, U.x((— °°,0)) = 
0, more explicit expressions are available for K and K. In this situation we take L t = 
— info<i<fXj, so the inverse process L7 1 is the passage time subordinator 

T v :=inf{f>0: inf X s < -y\, y > 0, (18) 

0<s<t 

and H, = t on {% < °o}. Let 

4>x(5)=inf{/3: yx(j3) = 5}, 5 > 0. (19) 

Since y/x is strictly decreasing on (— °°,<t>x(0)], the function <I>x : [0,°°) — > (— oo,<5x(0)] is 
the inverse of the restriction of y/x to (— °°,<I>x(0)]. It follows that 

K X (8,P) = p-3>x(8), 8>0,peR, (20) 

and as a consequence of the choice of normalisation in (17), 

* ) = |W' *>0,£>0. (2D 

See Section 8.1 of [25] or Section VII. 1 of [4] which, note, both apply to spectrally negative 
processes. If in addition g ,5f^ a \ then (21) remains true for 5 > 0, j3 > —a. 



3. Properties ofB and consequences for the insurance risk process 

We now present some analytical properties of the function B and discuss their implica- 
tions for the Levy insurance risk process. 
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m= , s ::r^" 2 : „v (23) 



3. 1 . Laplace transform and rate of growth of B. Direct analytic evaluation of B through 
(11) is not feasible so we turn to evaluating it by numerically inverting its Laplace trans- 
form. The following proposition provides the required theoretical result. 

Proposition 1. Assume (5) holds. Then for 8 > y/x(a) VO, 

B(8) := r e - St B(t)dt = ^ ~7~Tr~Ts \ • ( 22 ) 

Jo 8(8- v / x(«))K'(5,-a) 

When, further, X is spectrally positive, this takes the form 

<S>x{8)-a 

(8- Vx (a)) 2 <P x (8) 

Proof. Substituting for B from (1 1) we obtain 

8 r e - s, B(t)dt^8 fY 5 ' [' e^W'EeVX'-'dsdt 
Jo Jo Jo 

= 8 r ' e - s 'e^ a )' f e-Vx(^Ee a ^dsdt 
Jo Jo 

= S f°° J" ' e-l s -**W) t dte- 1 l*W'Ee cig 'ds (24) 

= T 1 —— [°° Ee al[ °8e- 5s 6s 

8 - Wx(a) Jo 

1 Ee aW ° 



where e is distributed as exponential with parameter 8 independently of X. Thus (22) 
follows from (16), and then (23) from (21). □ 

Next we investigate how (6) and complementary conditions relate to the growth of B. 
This will have potential consequences for modelling the insurance risk process. 

Proposition 2. Assume (5) holds, 
(i) IfEe aX ^ < 1 then 



(ii) IfEe aX * > 1 then 



S H = T~^~17\ \ e (°>°°)- (25) 

]im^^- = Yx(a). (26) 
(Hi) IfEe aX[ = 1 then (26) continues to hold (in which \j/x((x) = 0), together with 

liminf^ > 1; (27) 

thus, B(°°) = °°. If in addition EX\e aXi < °°, then 

B(t) 

limsup^-^ < °°. (28) 

t^oo t l 

Proof. Assume Ee aXl < 1. In that case yx((*) < 0, K (0,0) = q > and fc(0, — Of) > 0; see 
Proposition 5.1 of [24]. First integrate by parts and then use (22) to obtain 



e- dt B(dt)= 8e- 6t B(t)dt=- ^L> -, S>0. 

' (8-\j/ x {a))K(8,-a) 



Letting 8 1 then gives (25). 
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Now assume Ee aX ' > 1 . Observe that lnEe aXs is subadditive, hence by Fekete's lemma, 

]nEe aX ~* 

lim = r (29) 

for some r < °o. Since lnEe aXs > lnEe aXs — s\ir x {a), it follows that r e (a),°°). If 
r > \\fx(a), choose 8 e (\jfx(cc),r). It follows easily from (11) that 

liminf ^ > 8. (30) 

But then B(5) = °° which contradicts (22), so r = ty/x(a). 

Finally assume Ee 1 = 1. Then yx(oc) = so (27) is immediate from (11). In general 
whennj 6 J^"', we know only that Ee aX[ < °°, but now assume further that E(X\e aXl ) < 
°°. Since e aX ' is a (sub)martingale, by Doob's L 1 -maximal inequality (see Exercise 5.4.6 
of [14]), 

Ee aX " < (l-e- l ){l+E{aX+e aX »)). (31) 

Now 

X+e aX " < £(AXi) + n] =1 e aAX J (32) 

where AX, = Xj — i . Thus 

£(X+e aZ ") < n£ ((AXO+e"^ 1 ) {Ee a ^)"- 1 . (33) 

Hence substituting into (31) and using monotonicity, for some constant C < °° and all s > 

Ee aX s < c(1 + ^ e Vx(«)-v = +s ) ; (34) 

from which (28) follows after substitution in (11). □ 

Remark 1. From (12), the limit in (8) may be alternatively expressed as in (25). This 
is the form of the limit given in [24]; see Theorem 4.1 and Proposition 5.3 therein. The 
assumption EX\e aXx < °o arises in connection with Cramer's large deviation estimate for 
the probability of eventual ruin in the Vx((Z) = case; see [5]. 

If Ee aXl > 1 then (26) shows that B{t) grows exponentially with t . Thus for appropri- 
ately large u and t , increasing the time horizon by one unit increases the probability of ruin 
by a factor of e r , where r is essentially yfx(a). For example if r = 3 this is a factor of at 
least 20. This is clearly a situation which would concern any insurance company. 

If, instead, Ee 1 - 1 we are in the realm of the classical Cramer condition where B 
grows subexponentially. With the additional mild assumption EX\e aX ' < oo, B grows at 
most quadratically, as shown by (28). However B is still unbounded and the estimate in (8) 
is not uniform over all t > in this case. Further, the probability of eventual ruin is of a 
different order to the probability of ruin in finite time. 

If Ee aXl < 1 then none of the above issues arise. B is bounded, the estimate is uni- 
form in t and the finite ruin time probabilities are comparable to the infinite horizon ruin 
probabilities. Furthermore the modified form (13) of the limit holds, which as we will 
demonstrate provides a superior estimate for small u. 

In conclusion, unless exponential growth of the finite horizon ruin probabilities is to be 
modelled, then it is necessary that the claims surplus process satisfy Ee ' 1 < 1. Within 
this class, it may be desirable to further restrict to Ee aXl < 1 due to the intuitive appeal, 
and uniformity in f, of the estimate in (13). 
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3.2. Moments and smoothness of B. In this section we give some subsidiary results 
which expand on the properties of B. Since it will not be used in the remainder of the 
paper in an essential way, the proof of Proposition 3 is deferred to an Appendix. 

Proposition 3. Assume X is spectrally positive and has no Brownian component, (5) holds, 
and Ee aXl < 1 (so that\m\ t ^,ooX t = —°° a.s.). Then J^tB(dt) < °°; thus, the limit distribu- 
tion corresponding to B has finite expectation. 

To conclude this section, we mention some smoothness properties of B. Rewrite (1 1) as 

Bit) = e Vx W [' e-**( a) 'Ee cS 'ds. (35) 
Jo 

From this we see that B has a density B' satisfying 

B'{t) = \j/ x (a)B(t)+Ee aJ ', t >0. (36) 

If y/x (a) = 0, it follows immediately that B'(t) -> 1 as t -> and so B(t) increases approx- 
imately linearly near 0. The same conclusion holds when \j/ x (a)^ 0. For this first observe 
that from (35), 

L <B(t)<- —±Ee aX >. (37) 



Hence by (36), 



yx{a) yx(a) 



e Wx(a)t _ j + Ee aX, e yx(a)t Ee aX, 

< — )Ar < r~, — • (38) 



Vx(a) Vx(a) Vx(a) 

Thus again B'(t) tends to 1 as t — > 0. 



4. Tempered Stable Processes 

In this section we set out a parametric class of tempered stable Levy processes which 
will be used as the basis for later calculation and simulations. All processes will be spec- 
trally positive, i.e., have no downward jumps, but initially are not required to drift to — «>, 
so we denote them by Y to distinguish them from insurance risk processes which we will 
continue to denote by X. X will be obtained from Y in Section 5. 1 by subtracting a drift. 

4.1. Tempered stable processes. Let Y be a Levy process with characteristic triplet (^a^riy), 
where 

rr=c(l- P )- 1 -cj^ (l-e-™)^, a F 2 = 0, n F (dx) = — ^ *>0. (39) 

Then Y is a tempered stable process with parameters c > 0, a > and p e (0, 1) U (1,2). 
The choice of jy is made so that the cumulant of Y is 

\j/ Y (9) =lnEe eYl = -cT(-p)[a p - (a - 9) p ], 9<a, (40) 

where T denotes the usual gamma function. From this we find 

EYi = v4(0) = -cpT(-p)a p ~ l . (41) 

Y is a pure jump subordinator when p g (0,1) while for p <G ( 1 , 2) it is spectrally positive 
but of unbounded variation. 
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4.2. Inverse Gaussian processes. Suppose Y has characteristics given by (39) specialized 
by taking p = 1/2. Then Y is an inverse Gaussian process. It is a pure jump subordinator 
with cumulant 



and mean 



yy(0) = 2c^[%{s[a - yja - 0), < a, (42) 

EYi=cJ-. (43) 
V a 

5. Tempered Stable Insurance Risk Model 

We turn now to insurance risk modelling based on a Levy process X, obtained from 
the tempered stable process Y of Section 4.1, by subtracting a drift. We will focus on 
the case p € (0,1). We continue to require X t — > — °° a.s. to reflect that the insurance 
company intends to collect sufficient premiums to avoid certain ruin. A sufficient condition 
that ensures this is Ee aXl < 1. We note however that X t — > — °o a.s. may still hold when 
Ee aXl > 1 . This has some interesting consequences for insurance model parametrisation 
and safety loading management. 

5.1. Aggregate claims and the claims surplus model. Let Y be the tempered stable pro- 
cess of Section 4.1 with p e (0, 1). In that case Y is a pure jump subordinator which we 
can take to model the aggregate claims process. We may then consider a one parameter 
family of claims surplus processes indexed by the premium rate p: 

X t {p) =Y t -pt. (44) 

We will use a superscript p to denote quantities computed from X^ p \ for example its cu- 
mulant (0) — yfy(0)— p9. By the strong law for Levy processes, 

X t {p) -> -oo a.s. p > EY\. (45) 

Note also that by (40) and (41) 

(p) EY 

Ee aX i <1 ^=> \j/ Y (a)<pa <^ p > -cY{-p)a p - x p> — -, (46) 

P 

thus confirming that Ee aX i < 1 implies x} p ^ — > — °° a.s., via (45), since p € (0, 1). 

Taking p = 1 /2 we obtain a general model where aggregate claims are modelled by an 
inverse Gaussian process. The use of the inverse Gaussian process in this context is dis- 
cussed in Garrido and Morales [17], Morales [29] and Chaubey, Garrido, and Trudeau [9] 
among others. The choice p = 1/2 makes a number of computations easier. For example, 

it is a tedious but simple matter to confirm that $>x , defined in (19), is given by 



2nc 2 + 20TC ( yj {y/ap - y^Kc) 2 + 8p - y/(Xp ) + 8p 



-P 2 



5.2. Safety loading. In the Cramer-Lundberg model (2), if we write 

p = (l+£)AM, (47) 

then 4 is called the safety loading and, in practice, its value is typically of order 0.2 [18]. 

For the tempered stable model of Section 5.1, the natural interpretation of the safety 
loading is to write 

p=(l + $)EY 1 . (48) 
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Thus by (46), 

v4 p) (a)<0 P>Y^ ^ >l ~p~- (49) 

It is interesting to note that this condition imposes no restrictions on the parameters c and 
a. This will be further elaborated on in the next section. 

It has been suggested that the inverse Gaussian process be used as a model for aggregate 
claims. However when p = 1/2, % = 0.2 and p is given by (48), we are forced by (49) 
to consider the situation in which yA^(oc) > 0. Indeed this will be the case for any safety 

loading % < 1. Thus to obtain a model with a realistic safety loading and have y/^ (a) < 0, 
which, from Section 3.1 is necessary to prevent long term exponential growth of the finite 
time ruin probabilities, the inverse Gaussian process cannot be used. Instead, from (49), for 
any % > 0, to ensure that y/^ (a) < 0, where p is given by (48), one should take a tempered 
stable process with p e ([1 + 1). We will thus focus our numerical investigation of 
the ruin time estimates in (10) and (13) on processes satisfying this condition. 

This potentially undesirable aspect of the inverse Gaussian process results from the as- 
ymptotic (in t) behaviour of the asymptotic (in u) estimate (11), together with the above 
safety loading considerations. For small values of the initial reserve or over short time pe- 
riods, exponential growth in the finite time ruin probabilities may not be exhibited. In this 
case the inverse Gaussian process with a safety loading of 0.2 may prove to be an adequate 
model. Note however that estimate (13) is not available in this case since ^ x p \a) > 0. 

5.3. Interpretation of model parameters. We next discuss the role of the various param- 
eters in these models. Fix a,b>0 and let R t = bY at . Then R represents the same aggregate 
claims process but with different units of currency and time. For example if Y, is the aggre- 
gate claims after t years measured in millions of dollars, and a = 1/4 and b = 10 3 , then R t 
is the aggregate claims after t quarters measure in thousands of dollars. The Levy measure 
ofR is 

„ , , , abPce-^l^dx 

U R (dx) = — p , x > 0. (50) 

This is of the same form as Yly but with different values for the parameters c and a. Thus 
varying c and a is equivalent to changing the currency and time scale. This is not the case 
for p , though. Similarly if we write p = ( 1 + 1, )EY X , so that X t (p) = Y t - ( 1 + E, )EY t , then 

bX^=R t -{\ + S,)ER u (51) 

and we see that the safety loading also does not depend on the units of currency or time. 
Thus, up to a change of scale, the key parameters to vary to obtain different models are p 
and |. 

6. Numerical Approximation 

This section introduces the techniques to be used in numerically approximating the ruin 
time distribution, when X is spectrally positive, using (10) and (13). 

6.1. Approximating B. Expressions for the Laplace transform of B are given in Proposi- 
tion 1, from which B(t) can in principle be evaluated using the Bromwich integral 

B(t) = — e 5t B{8)d8, (52) 

27T1 J £ -i<x, 
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where £ is chosen so that B in (22) is analytic in the region 91(5) > £. An extensive 
coverage of methods for approximating integrals of this form is presented by Cohen in 
[10]; in particular, benchmarks comparing the relative errors and computational times for 
a number of the methods are in [10]. As a check, we used two different approaches for the 
calculation of (52). 

Fixed-Talbot method. The first consists of a straightforward approach by Valko and Abate 
[34] which relies on multi-precision arithmetic called the fixed-Talbot method; see [10], 
p. 138. This method is very attractive from an implementation point of view if one has 
access to a software package with arbitrary precision arithmetic capabilities (e.g., Mathe- 
matica [28]). The fixed-Talbot approach consists of deforming the contour in (52) to the 
path 8(9) = r0(cot(0) +i) for — n < 9 < n, where r is a parameter. Integration over this 
new contour gives 

B(t) = — I' e tS ^B(8(9))8'(9)d9. 

Substituting S'(0) = ir(l +iff(0)) with a{9) := 9 + (0cot(0) - l)cot(0) and knowing 
that B is real-valued, we obtain 

B{t) = T - J* 3i(e t5 WB(8{9))(l +itr(0))) d9. 

This integral is then approximated using a trapezoidal rule with step size n/M and 9j = 
jn/M by 

B M (t) = - ( \B{r)e rt + - £ 9l((l +icr(0 J )) J B(5(0 J )) e 5 ^')J . 

As suggested by Valko and Abate, we choose r = 2M/(5?) where M is the number of 
decimal digits of precision required. 

Levin method. The second method provides an alternative for situations where arbitrary 
precision arithmetic is not available. It starts from the observation that B(t) = for t < 
and B is real-valued, so (52) simplifies to 

2e £t r ~ 

B(t) = / 3i(B(e + iu))cos(ut)du. (53) 

7T Jo 

This is an integral of oscillatory type so care has to be taken in performing a numerical 
approximation. We specialised the approach of Levin [27] to our particular case (53), as 
follows. 

For t > and a large enough M > we can approximate (53) by 

2 f M 

B(t) w — / f(u) cos(tu) du 
% Jo 

where f(u) := 3iB(iu). Assume /(«) to be of the form 

f(u)=F'(u)-tF(u)tm(tu), 0<u<M, (54) 
for some function F. Then it follows that 

i-M i-M 

I f(u)cos(tu)du — / (F'(u)cos(tu) — tF(u) sin(fw)) du 
Jo Jo 

= / — (F(u)cos(tu)) du 
Jo du 

= F(M)cos(tM)-F(0). 
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FIGURE 1 . B (8) from (22) and the corresponding B(t) for t > 0, by nu- 
merical inverse Laplace transform. Generated from the tempered stable 
process in (44) with Y having triplet (39), where p = 0.99, a = 1 .0, 
c = 0.01, andt, =0.2. 

To find this unknown function F, we consider relation (54) as an ordinary differential 
equation (ODE) where / is given and F is to be determined. Since a solution of (54) seems 
difficult to obtain in closed form, we solve the ODE using a numerical method. Assume 
there exists n EN large enough so that F can be approximated arbitrarily closely on the 
interval [Q,M] by a function F„ given by a linear combination of n predetermined basis 
functions pi(u),...,p n (u): 

n 

F n(u)=Y, c kPk(u), (56) 
k=\ 

where c\, . . . ,c„ are n unknown coefficients to be determined. For simplicity, we made 
the choice pk(u) = 71(m), where T^iu) is the £-th Chebyshev polynomial of the first kind 
(i.e., Tk satisfies 7^(cos(w)) = cos(fcw)). Substituting (56) into (54) and using the identity 
TL{u) = kUk-\{u), where U^u) is the Chebyshev polynomial of the second kind, we obtain 
the following equation 

n 

f(u) = £ c k (kU k -i («) - tT k {u) tan(f«)) . (57) 

k=\ 

To find the coefficients ci,...,c„, we choose n collocation nodes = U\ < «2 < """ < 
u„ = M and evaluate (57) at these points in order to set up a system of n equations with n 
unknowns which can be solved for the c,. Once these coefficients are obtained, and since 
7yt(0) = cos(^7r/2), we have the approximation 

B ( f ) ~ | ^LQ(7i(M)cos(fM)-cos(^/2))^ . (58) 

Here we see that the approximation depends on good choices of the cut-off value M > 0, 
the number n of basis functions used to approximate F, and the location of the colloca- 
tion nodes u\,...,u n . After inspection of the behaviour of B near zero (see for example 
Figure 1), we chose M = n and set the location of the nodes at m,- = cot(/7r/(2n)), so they 
accumulate near zero. 
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6.2. Approximating Tl x (u) and P{l{u) < °°). Approximating Tl x (u) for an arbitrary 
choice of u > is straightforward. For example, if X is a tempered stable process then 

— + r°° ce~ ax 

nj(«)= / -^fc, ">0 (59) 

(see (39)). The integral here is easily calculated using a numerical quadrature routine, or a 
simple asymptotic approximation may be adequate for a large enough u: 

— + ce~ au 

n ^ (M) = ^T^ (1+ ° (1)) aSM ^~' (60) 

For the models we consider, II X (m) has a singularity at 0, so it could be expected that 
the estimate for P(x(u) < t) given by (10) is large for small u. This suggests that (13) 
may provide better estimates, at least for small u, when applicable. When X is spectrally 
positive, it is well-known that the infinite horizon ruin probabilities P(x(u) < °°) can be 
approximated arbitrarily closely by numerically inverting a Laplace transform; specifically, 
if EX\ < 0, then from Theorem 8.1 of [25] we have that 

P(t(m) < oo) = 1 +EX l W(u) (61) 

where the scale function W satisfies 

/ e-P u W(u)du = (62) 

Jo Wx(-P) 

Recently, an extensive overview of scale functions and their numerical evaluation through 
the inversion of a Laplace transform has been presented by Kuznetsov, Kyprianou, and 
Rivero [26]. 



7. Simulation Methodology 

Simulation of observations on random variables with tempered stable distributions and 
the resulting processes has recently become an active topic of research due to their use in 
a variety of different applications; see the recent survey by Kawai and Masuda [23] and 
the references therein. However, to the best of our knowledge, simulation of finite-time 
ruin probabilities for tempered stable processes has not yet been covered in the literature. 
Consequently, we provide some detail on our approaches. 

Naive approach. A naive approach is tempting; simulate the tempered stable process in- 
crement by increment by sampling from the increment distribution and tally the number of 
paths that pass above level u by time t, then calculate the ratio of crossing paths to total 
paths simulated. When the stability index p < 1, the law of the process increments can be 
simulated exactly. On the other hand, when p > 1 no practical exact simulation method 
exists and one must resort to an approximation [23]. This is particularly troublesome when 
the probability of ruin is very small due to the possibility of bias in the convergence of 
the simulation algorithm. In addition, as the simulation of tempered stable random vari- 
ables is currently not built-in to standard software packages, a custom implementation is 
needed. Further, simulation of a tempered stable random variable attracts a computational 
cost greater than that of a stable random variable. These disadvantages motivated us to 
develop another approach that we shall now explain. 
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Measure change approach. There is a useful relationship between tempered stable and 
stable processes that we can exploit. Consider a spectrally positive stable process Z = {Z t : 
t > 0} of index p e (0, 1) U (1,2) with characteristic triplet (y z ,0,IT z ), where 

c c dx 

72 = 7— -, n z (ck) = - rT ^, x>o,oo. (63) 

It is easily shown, by integration by parts, that for p e (0, 1) U (1,2), X > 0, 

/ («-**- 1 +AxL {0< . v<1} )n z (dx) = cAT(-p) + — -. (64) 
Jo 1 ' 1 - p 

Hence the Laplace exponent of Z\ satisfies 

y/ z (-A) =log£e _AZl =cr(-p)A p , A>0. 

When p e (0, 1), the process Z is a pure jump subordinator while for p <E (1,2) it is spec- 
trally positive with finite mean but of unbounded variation. 

Let z\ p) = Z, pf . The characteristics of and Z (p) are - p,0,IIy) and (y z - 
p,0,n z ) respectively. Now may be obtained from by an exponential change of 
measure, i.e. is an Esscher transform of X^ p \ Specifically assume XW and zW are 
given on a filtered probability space (£2,^, (^ t )t>o,P) and define 

g= e ^ W -^W onJ r f . (65) 

Then under Q has the same characteristics as Z^ ; see Theorem 2 of VII. 3c in [32]. 
Rewriting (65) as 

dP = e- ca F >+1 # ) WdQ on& t (66) 

we find 

P{x x(p) («) < f) = E Q (e-«*P + &W; t x(p) (u) < t) 

= E(e- aZ > P) ;T z{p) (u) < r) e V* W (67) 
= E(e- aZ > P) ; T zip) („) < t) e -W-P)«"+P«> 

Thus to calculate ruin probabilities for X^\ we need only simulate the stable process zM. 
This has a number of advantages. First, numerical packages that simulate stable random 
variables are readily available. Second, the simulation of an increment of a stable process 
is less computationally expensive than simulating a tempered stable process. And third, 
the law can be exactly sampled in the case p > 1 . 

Implementation details. We construct the stable process ZW from a samples of a stable 
random variable, S, as follows. Suppose available a numerical package for simulating a 
general stable random variable S ~ Stable (p , j3 , ji , V ) with index p € (0,1)U(1,2), skew- 
ness parameter p\ location parameter /x, and scale parameter v, having characteristic expo- 
nent expressed in the form (e.g., [31, Equation (14.24)]) 

^(0) = -v|0| p (l -ij3sgn0tan^) +i^0. (68) 

Then (cf., e.g., proof of [31, Theorem 14.10]) the law of Z^' 5 ', for h > 0, is obtained from 
the law of S by setting 

jS = 1, jti. = -ph, v = -chcos(7cp/2)r(-p). (69) 
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We simulate a path t M> Z, = Z t — pt by decomposing the path as a sum of increments 
over small time intervals h > 0; see Algorithm 1. 

Algorithm 1 P(t(u) < t ) using time increment h > and n simulations 

1: sum <— 
2: H < ph 

3: v<- (-/iccos(7rp/2)r(-p)) l /P 
4: for / = 1 ,...,« do 

5: S <- 0, X <- 

6: hit ■<— false 

7: while 5 < f do 

8: dX~StabIe(p,l } |U,v) 

9: S •<— 5 + ft 

10: X^X + dX 

11: ifX>«then 
12: hit «- true 

13: if hit then 

14: sum = sum + e~ aX 

15: return (sum/n)exp(— (cT(-p)a p +pa)t) 



To get a good approximation of the ruin probability for a given (u,t) when P(t(u) < t) 
is small, a large number of simulations n and a small time step h > are needed. To speed 
up the simulations, we implemented Algorithm 1 in C++ 11 and modified the algorithm 
slightly to run in parallel using DpenMP. In the parallelized version of our algorithm, each 
thread received its own copy of a random number generator (mtl9937: Mersenne twister) 
initialised with independent initial seed. Since the stable distribution is not part of C++ 11 
we wrote our own implementation based on an acceptance-rejection method that requires 
only uniform variates and exponential variates (both available in C++11), Zolotarev's func- 
tion, and the function sinc(jt) := sin(x)/x [11]. This implementation takes the random 
number generator as an argument so that the generation of variates is independent across 
threads. The for loop at line 4 in Algorithm 1 is distributed over the threads using an 
OpenMP parallel pragma and the shared variable sum is set to be reduced using the addi- 
tion operator. 

In Table 1 and Table 2, we simulated P(t(u) < f) at u =0.1 and t = 2.0 using the 
naive and the measure changed algorithm with step sizes of h = 0.0001 and h = 0.01, 
respectively. We also fixed p = 0.99, c = 0.01, a — 1 and | = 0.2. The tables contain the 
time 3 in seconds taken to simulate P{t(u) < t) using n sample paths, the mean value of 
P(x{u) < t) calculated using N — 30 batches of n sample paths, and a/y/N, where c is the 
standard deviation of P(i(u) < t) over the batches. 

The measure change algorithm performs better than the naive algorithm as it results in a 
smaller confidence interval (i.e., a/y/N is smaller) and shorter run times for the parameter 
values we are interested in. In some cases, the naive approach can give shorter run times. 4 
The main advantages of the measure change approach are the ability to accurately simulate 
the case p > 1 and the simplicity of implementation as it does not require the simulation 
of tempered stable random variables. 



3 On an Intel Xeon W3680 @ 3.33Ghz using 12 threads. 

4 This can occur, for example, when t is large, as the naive approach does not need to simulate the full path up 
to time t if the process passes above u at some earlier time. 
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Naive 






Measure Chan 


? e 


n 


Time (s) 


Mean 




Time (s) 


Mean 




32 


0.079217 


0.051041667 


0.0085395638 


0.076420 


0.057849267 


0.0064226185 


64 


0.153592 


0.0484375 


0.0056237028 


0.141648 


0.0543778 


0.0038164514 


128 


0.280040 


0.054947867 


0.0048331632 


0.264577 


0.052818533 


0.0031319083 


256 


0.558673 


0.051822867 


0.0027775502 


0.526003 


0.052552067 


0.0022664163 


512 


1.080112 


0.052604233 


0.0016828511 


1.040853 


0.052797033 


0.001440741 


1024 


2.173619 


0.051692733 


0.0011881599 


2.240990 


0.051077267 


0.00091771839 


2048 


4.281521 


0.049755767 


0.0010048818 


4.147853 


0.050568033 


0.00064971641 


4096 


8.624235 


0.0508301 


0.00068509617 


8.315087 


0.050364233 


0.00050435188 


8192 


17.148903 


0.0512859 


0.00049187434 


16.571840 


0.050474467 


0.00036946037 


16384 


34.222089 


0.0512817 


0.00032080185 


33.201601 


0.050540433 


0.00029116668 



TABLE 1. Mean of P(x{u) < t) calculated using the two approaches, 
calculated using N — 30 batches ofn paths and a time step ofh = 0.0001. 







Naive 






Measure Chan 


ge 


n 


Time (s) 


Mean 




Time (s) 


Mean 




32 


0.004348 


0.046875 


0.0069877124 


0.004592 


0.046198033 


0.0057320205 


64 


0.005102 


0.044270833 


0.0057095406 


0.004866 


0.046170833 


0.0039138265 


128 


0.007268 


0.045833433 


0.0039947614 


0.006323 


0.047625433 


0.0028581633 


256 


0.009335 


0.046354167 


0.0021084348 


0.009371 


0.047685367 


0.0017811631 


512 


0.015461 


0.049869767 


0.0014499284 


0.014557 


0.0495607 


0.0013563884 


1024 


0.029704 


0.0511068 


0.0011481938 


0.024450 


0.0504672 


0.00093407293 


2048 


0.047680 


0.050992767 


0.00093047734 


0.044822 


0.050154467 


0.00067853448 


4096 


0.090405 


0.050187133 


0.00061191755 


0.093858 


0.050157833 


0.00042720285 


8192 


0.178851 


0.050354 


0.00035644764 


0.167700 


0.050329367 


0.0003197472 


16384 


0.344298 


0.050529 


0.00026127832 


0.329710 


0.050477933 


0.00029453508 



TABLE 2. Mean of P{x{u) < t) calculated using the two approaches, 
calculated using N = 30 batches ofn paths and a time step ofh = 0.01. 



8. Comparison of Asymptotic and Simulation Estimates 

In this section, we report on the estimates of the ruin time distribution obtained from 
the asymptotic formula (10), and from the modified version (13), and compare them to the 
values obtained from the Monte Carlo simulation as described in Section 7. 

As shown in Figure 2, it may happen that the leading term in formula (10) overestimates 
the value of P{l{u) < t) for small u. Direct application of the asymptotic estimate (10) can 
give putative values of P(t(u) < t) greater than 1 when u is small. This, of course, is a 
result of the singularity in TY^{u) at u = 0. However, this is not an issue with the modified 
estimate (13), which provides meaningful estimates even for very small u. 

The value of in (13) can be calculated by combining (25) and (21), together with 
the observation that q = \EX^\. This yields 

a\Ex\ p) \ 



B i°°) = _lpp__L. (70 ) 

Evaluation of P{x{u) < °°) is through (61) and (62) using the algorithms for numerically 
inverting Laplace transforms discussed in Section 6. 1; see Figure 3. Substituting the result- 
ing values into (13) leads to a substantially improved estimate. To benchmark the estimate 
(13), we compared it against the Monte Carlo simulation. Parameter values were chosen so 
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FIGURE 2. The value of Yl x (u)B(t) for various (u,t) [left] and for 
t = 1,5, 10 [right]. Note that values greater than 1 occur. Model and 
parameter values as in Figure 1. 




FIGURE3. The scale function ut-^W (u) in (62) [left] and u i— ► P(z(u) < 
°°) in (61) [right]. Model and parameter values as in Figure 1. 



that the resulting model reasonably well approximates practice. In this regard, we followed 
Grandell [19] p. 145 who cites working actuaries as regarding a safety loading of 0.2, an 
initial reserve equal to the expected aggregate claims over a one year time period, and a 
planning horizon of five years, to be practical. For added flexibility, we allowed for an 
increased planning horizon of up ten years. 

As discussed in Section 5.3, the key parameters in the tempered stable models are p 
and 4 . The parameters c and a in (39) correspond to changes of scale. In the scenario 
illustrated in Figure 4, we took a safety loading of % = 0.2. In order that (49) hold we then 
must choose p £ (5/6, 1). After some experimentation we found the asymptotic estimate 
performs better for values of p close to 1 . We took p = .99. The values of c and a may then 
be chosen to set a convenient scale in t and u. For example in Figure 4 we took c = .01 and 
a = 1. We assume one time unit corresponds to 6 months. Then the expected aggregate 
claims over one year is EY2 = 1.9886. Thus we plot the ruin probabilities for < t < 20 
and < u < 2. 5 5 



To illustrate the use of c and a in setting the scale, if we took c = (.01)2 °' and a = 2 then the plots in 
Figure 4 would change only in the scale on the horizontal axes, which would become < t < 10 and < u < 1. 
From Section 5.3 this corresponds to ruin probabilities for the new process R t = X% /2 where is the original 
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FIGURE 4. P(t(w) < t) for Monte Carlo (magenta) and asymptotic for- 
mula (cyan) [left] and the ratio of the two [right]. Model and parameter 
values as in Figure 1. 



As shown in Figure 4, the estimate (13) performs very well in this case when t is greater 
than 10 (corresponding to 5 years) and u is greater than 1 (corresponding to half the ex- 
pected annual aggregate claims). As Table 3 indicates, when u = 2 and t > 10 the relative 
errors are less than 8 .2% and decrease to below 2% by the time t reaches 20. Even for u = 1 
the relative error is less than 14% for t > 10. The variability in the simulated probabilities 
observed in Table 3 arises because the probability being estimated is very small. This 
accords with an insurance company's desire to set its ruin probability over the planning 
period to be negligible. In this example it is of order 10~ 3 which seems reasonable. 

In Table 3, the infinite horizon probabilities also provide reasonably good estimates for 
the simulated probabilities. In general, P(l(u) < °°) can always be used as an upper bound 
for the probability of ruin in finite time, but there is a question of how precise a bound 
it may provide. It is quite possible that it may grossly overestimate the finite time ruin 
probability. We observed empirically, for a wide range of parameter values, that (13) gives 
a lower bound for the probability of ruin in finite time calculated via simulation, and when 
it does not, it overestimates only slightly. Thus (13) combined with the infinite horizon 
ruin probability can be used to place good bounds on the finite time ruin probabilities. 

9. Conclusion 

Up till the publication of the estimate (10) in [21], simulation has been the only method 
of calculating the distribution of the ruin time in the convolution equivalent model. Our 
aim in the present paper was to show that the function B(t ) can be calculated numerically 
in an interesting and useful class of models, and to examine some of its properties. The 
formula (10), though asymptotic, is fast to calculate, and is immediately useful for initial 
calibration of a model. Once initial estimates of parameters are found in this way, a full 
Monte Carlo simulation could be performed to refine the estimates if desired. Additionally, 
the function B(t) in (11) or its normalised counterpart in (13) can be analysed to provide 
much insight into the ruin time distribution in these models. In particular it provides new 
insight into safely loading management for these models. Finally we observed empirically 
that (13) seems to give a lower bound on the probability of ruin in finite time calculated via 



process in Figure 4 with c = .01 and a = 1. In this case t would now be measured in years, and the expected 
annual claims would be £7?i = EY 2 /2 = 0.9943. 
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u 




a 


s 


i 


a/s 


i/s 


\a — s\/s 


| i — s\ / s 


1 


10 


0.00330802 


0.003835 


0.00393118 


0.862586 


1.02508 


0.137414 


0.0250804 


1 


12 


0.00350162 


0.003829 


0.00393118 


0.914499 


1.02669 


0.0855013 


0.0266867 


1 


14 


0.00363522 


0.003872 


0.00393118 


0.938849 


1.01528 


0.0061151 


0.0152849 


1 


16 


0.00372736 


0.003725 


0.00393118 


1.00063 


1.05535 


0.000634196 


0.0553512 


1 


18 


0.00379087 


0.003704 


0.00393118 


1.02345 


1.06133 


0.0234522 


0.0613346 


1 


20 


0.00383461 


0.003971 


0.00393118 


0.965654 


0.989973 


0.0343456 


0.0100269 


1.5 


10 


0.00127612 


0.0015 


0.00151651 


0.850746 


1.01101 


0.149254 


0.0110093 


1.5 


12 


0.0013508 


0.001506 


0.00151651 


0.896947 


1.00698 


0.103053 


0.00698137 


1.5 


14 


0.00140234 


0.001573 


0.00151651 


0.891509 


0.96409 


0.108491 


0.0359098 


1.5 


16 


0.00143789 


0.00138 


0.00151651 


1.04195 


1.09892 


0.041947 


0.0989231 


1.5 


18 


0.00146238 


0.001549 


0.00151651 


0.944083 


0.979028 


0.0559169 


0.0209723 


1.5 


20 


0.00147926 


0.001395 


0.00151651 


1.0604 


1.08711 


0.0604019 


0.0871068 


2 


10 


0.000543995 


0.000585 


0.000646473 


0.929906 


1.10508 


0.0700941 


0.105081 


2 


12 


0.000575831 


0.000627 


0.000646473 


0.918391 


1.03106 


0.0816087 


0.0310568 


2 


14 


0.000597803 


0.000648 


0.000646473 


0.922535 


0.997643 


0.0774648 


0.00235709 


2 


16 


0.000612955 


0.000574 


0.000646473 


1.06787 


1.12626 


0.0678655 


0.126259 


2 


18 


0.000623398 


0.000642 


0.000646473 


0.971025 


1.00697 


0.0289752 


0.00696668 


2 


20 


0.000630592 


0.00062 


0.000646473 


1.01708 


1.0427 


0.0170838 


0.0426978 



TABLE 3. Comparison and benchmark of asymptotic formula for finite- 
time ruin probability (a) against Monte Carlo simulation (s) and infinite- 
horizon ruin probability (i). Simulations were performed with Algo- 
rithm 1 using time increment h = 0.001 andn = 32768 trials. 



simulation, which combined with the infinite horizon probability P(t(u) < °°) as an upper 
bound, might have useful practical applications. 

10. Appendix 
Here we prove the proposition in Section 3.2. 

Proof of Proposition 3. Assume (5) and Ee aXl < 1. Let B(t) := B(°°) -B(t), t > 0, and 
use (22) and (25) to write 



/ e- Xt B{t)dt = X- l B{°°)- ( e- Xt B{t)dt 
Jo Jo 

q K (X,Q) 



(71) 



-\l/ x (a)K(0,-a) (X-Yx(a))K(X,-a) J ' 
The last expression can be simplified to 

-C l (X,a){K(X,0)-q)+XC2(X,a)+C 3 (X,a)(K(X,-a)-K(0,-a)) 

XC 4 {X,a) 

where 

limCi(A,a) =lim(A- y/x(a))?c(A, -a) = -\j/x(a)K(0,-a), 
\imC2(X,a) = limK(X ,0)k(X , -a) = qK(0,-a), 



]imC 3 (X,a) =]im-\i/ x (a)K(X,0) = -qyxia), 
x 4,0 x 40 



and 



limC4(A,a) = lim-y/x(a)K-(0,-a)(A - \j/x(cc))K(X,-a) = (-y/x(a)) k 2 (0,-a). 

x 4,0 x 40 
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All four limits are finite and strictly positive. Also note that 

C 3 (0,a)-Ci(0,a) = -qyx{a) + y/x(a)K-(0,-a) = ~Yx(a)(q- jc(0,-a)) > 0. 
Recalling (15), we get by monotone convergence 

r *(*.0)-g v Ad L - 1+ X >0 (l- e -^)n L -,(A) y 

hm — — t = hm =dr-i+/ fllr-i (at , (72) 

A|0 A A|0 A L Jr>o 

in the sense that both sides are finite or infinite together. Also observe that 

K(X,-a)-K(0,-a) = Xd,-i+ [ ( 1 -e~ lr ) U L -i (dt) 

7 ( >oV ) (?3) 

+ L( 1 -^)L( e "- 1 ) n - 1 .-^^- 

Dividing by A and letting A \. gives 

,. fc(A,-a) - K(0,-a) , /" _ , , . 

hm— — '— — — -=dr-i+/ tU r -i(dt) 

Xio A L Jt>o 



t I [e ah -\)Xl L -\„{dt,dh) 

t>0 Jh>0 



(74) 



again in the sense that both sides are finite or infinite together. Returning to (71), and 
letting A 4- we find that 



C 4 (0,a) / B(t)dt = (C 3 (0,a)-Ci(0,a)) ( d L -i + / rn L -,(^ 
Jo V ^'>o 



(75) 

■C 3 (0,a) / f/ (g a *-l)n £ -i ff (A,dA). 



r>0 V/i>0 



Since C 4 (0,a) > 0, C 3 (0,a) > Ci(0,a) and C 3 (0,a) > 0, we find that £ B{t)dt < °o if 
and only if the two integrals on the righthand side of (74) are finite. 

When in addition X is spectrally positive, the second integral is always finite, while if 
additionally Ox = the first integral is also. To see this, we first note that the integral 
f h> ie ah Tln(dh) is finite whenever Ee aXl < °o by Proposition 7.1 of [19]. Then, treating 
first the double integral in (74), we have for the component over < t < 1 , 



f tf (e ah -l)n L -i „{dt,dh) 
J0<r<l Jh>0 V ' 



< t [e a l {0<h<1} +e ah l {h>1} )n L -i„(dt) (76) 

J0<r<l Jh>0 ^ ' 

<e a ( tU L -i(dt)+ [ e ah Xl H (dh) < °°. 

J0«<1 Jh>l 

(Here the first integral on the righthand side is finite, as for any subordinator.) 

To deal with the remaining part of the integral over t > 1, we assume further at this stage 
thatX is spectrally positive. Thus by [25], p. 208, 

n L -i H (dt,dh) = I llx(dh + v)P(% e dt)dv, t>Q,h> 0. (77) 

Jv>0 
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Since X, — » — °° a.s. it follows from Theorem 1 of Doney and Mailer [13] (applied to — X) 
that Ex\ < oo. Thus 



t( (e ah -l)ll L - lH (dt,dh) 

Jh>0 V / 



t>l Jh>Q 

ah 



t [e an -\)l n x (dh + v)P(T v edt)dv 

r>l Jh>0 V / Jv>0 

,a(h-v) _ A Ux ( dh ) f t pfc g dt } dv 
v>0Jh>v\ ' Jt>l 

r r / \ 

</ / ( e a/ '-i) e -« i 'nx(^)£(T y )rfv 

Jv>0Jh>v V / 

= £?i / t ve- av dvU x (dh) 

Jh>0 V /JO 



<£Ti / li c""-l nj((i(i)+£T| ve- av dv e an H x {dh). 

J(Xh<\ \ ' JO Jh>l 

The first integral on the righthand side is obviously finite, and the second term on the 
righthand side is finite since Ee aXl < oo. 

Finally we deal with the first integral on the righthand side of (74). We need only 
consider values of t > 1. For this we further assume a x = 0. Since dg > it follows from 
Corollary 4 of [ 12] that d# = 0. Hence X does not creep up and so by Theorem 3.4 of [20], 
(77) also holds for h = 0. Thus 

n L -i (dt) = [ U x (v)P(% € dt)dv, t > 0, (79) 

Jv>0 

and hence 

fllr-i (dt) = / n^(v) / tP(% e dt)dv < ETi I vn^(v)c/y < °o 

?>1 Jv>0 Jt>l Jv>0 

since E(X^) 2 < °° as a result of Ee aX] < So the last integral converges too. 

Thus the two integrals on the righthand side of (74) are finite and so [ Q ™tB(dt) is finite. 

□ 
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